home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / DRAND.C < prev    next >
C/C++ Source or Header  |  1996-03-17  |  3KB  |  158 lines
  1. /* drand.c
  2.  *
  3.  * Pseudorandom number generator
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double y, drand();
  10.  *
  11.  * drand( &y );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Yields a random number 1.0 <= y < 2.0.
  18.  *
  19.  * The three-generator congruential algorithm by Brian
  20.  * Wichmann and David Hill (BYTE magazine, March, 1987,
  21.  * pp 127-8) is used. The period, given by them, is
  22.  * 6953607871644.
  23.  *
  24.  * Versions invoked by the different arithmetic compile
  25.  * time options DEC, IBMPC, and MIEEE, produce
  26.  * approximately the same sequences, differing only in the
  27.  * least significant bits of the numbers. The UNK option
  28.  * implements the algorithm as recommended in the BYTE
  29.  * article.  It may be used on all computers. However,
  30.  * the low order bits of a double precision number may
  31.  * not be adequately random, and may vary due to arithmetic
  32.  * implementation details on different computers.
  33.  *
  34.  * The other compile options generate an additional random
  35.  * integer that overwrites the low order bits of the double
  36.  * precision number.  This reduces the period by a factor of
  37.  * two but tends to overcome the problems mentioned.
  38.  *
  39.  */
  40. #include "mconf.h"
  41. #include <stdlib.h>
  42.  
  43. /*  Three-generator random number algorithm
  44.  * of Brian Wichmann and David Hill
  45.  * BYTE magazine, March, 1987 pp 127-8
  46.  *
  47.  * The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12.
  48.  */
  49.  
  50. static int sx = 1;
  51. static int sy = 10000;
  52. static int sz = 3000;
  53. static double unkans = 0.0;
  54.  
  55. /* This function implements the three
  56.  * congruential generators.
  57.  */
  58. void
  59. DrandSeed(unsigned UserSeed)
  60. {
  61.     sx = (UserSeed & RAND_MAX) | 1;
  62.     sy = (UserSeed + 10000) & RAND_MAX;
  63.     sz = (UserSeed +  3000) & RAND_MAX;
  64. }
  65. static    void
  66. ranwh(void)
  67. {
  68. int r, s;
  69.  
  70. /*  sx = sx * 171 mod 30269 */
  71. r = sx/177;
  72. s = sx - 177 * r;
  73. sx = 171 * s - 2 * r;
  74. if( sx < 0 )
  75.     sx += 30269;
  76.  
  77.  
  78. /* sy = sy * 172 mod 30307 */
  79. r = sy/176;
  80. s = sy - 176 * r;
  81. sy = 172 * s - 35 * r;
  82. if( sy < 0 )
  83.     sy += 30307;
  84.  
  85. /* sz = 170 * sz mod 30323 */
  86. r = sz/178;
  87. s = sz - 178 * r;
  88. sz = 170 * s - 63 * r;
  89. if( sz < 0 )
  90.     sz += 30323;
  91. /* The results are in static sx, sy, sz. */
  92. }
  93.  
  94. /*    drand.c
  95.  *
  96.  * Random double precision floating point number between 1 and 2.
  97.  *
  98.  * C callable:
  99.  *    drand( &x );
  100.  */
  101. void
  102. drand(double *a)
  103. {
  104. unsigned short r;
  105. unsigned short *p;
  106. #ifdef DEC
  107. unsigned short s, t;
  108. #endif
  109.  
  110. /* This algorithm of Wichmann and Hill computes a floating point
  111.  * result:
  112.  */
  113. ranwh();
  114. unkans = sx/30269.0  +  sy/30307.0  +  sz/30323.0;
  115. r = (unsigned short) unkans;
  116. unkans -= r;
  117. unkans += 1.0;
  118.  
  119. /* if UNK option, do nothing further.
  120.  * Otherwise, make a random 16 bit integer
  121.  * to overwrite the least significant word
  122.  * of unkans.
  123.  */
  124. #ifdef UNK
  125. /* do nothing */
  126. #else
  127. ranwh();
  128. r = (unsigned short) (sx * sy + sz);
  129. p = (unsigned short *)&unkans;
  130. #endif
  131.  
  132. #ifdef DEC
  133. /* To make the numbers as similar as possible
  134.  * in all arithmetics, the random integer has
  135.  * to be inserted 3 bits higher up in a DEC number.
  136.  * An alternative would be put it 3 bits lower down
  137.  * in all the other number types.
  138.  */
  139. s = *(p + 2);
  140. t = s & 07;    /* save these bits to put in at the bottom */
  141. s &= 0177770;
  142. s |= (r >> 13) & 07;
  143. *(p + 2) = s;
  144. t |= r << 3;
  145. *(p + 3) = t;
  146. #endif
  147.  
  148. #ifdef IBMPC
  149. *p = r;
  150. #endif
  151.  
  152. #ifdef MIEEE
  153. *(p + 3) = r;
  154. #endif
  155.  
  156. *a = unkans;
  157. }
  158.